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Characteristic Lyapunov vectors in chaotic time-delayed systems 
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We compute Lyapunov vectors (LVs) corresponding to the largest Lyapunov exponents in delay-differential 
equations with large time delay. We find that characteristic LVs, and backward (Gram-Schmidt) LVs, exhibit 
long-range correlations, identical to those already observed in dissipative extended systems. In addition we give 
numerical and theoretical support to the hypothesis that the main LV belongs, under a suitable transformation, 
to the universality class of the Kardar-Parisi-Zhang equation. These facts indicate that in the large delay limit 
(an important class of) delayed equations behave exactly as dissipative systems with spatiotemporal chaos. 
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I. INTRODUCTION 

Delayed dynamical systems (DDSs) serve to model diverse 
phenomena in physics lU (prominently in optics), but also 
in other fields 121 such as engineering, biology, climatology 
or ecology. More than one decade ago (3|-|5[] it was found 
that delayed systems with one constant delay can be studied 
like extended systems, and they present not only analogies 
but equivalent phenomena such as pattern instabilities |g] and 
high-dimensional chaos tM- Chaos in DDSs is considered as 
a particular type of spatio-temporal chaos for which the delay 
plays the role of the system size. 

There are a number of studies concerning the tangent dy- 
namics of systems with time delay. Previous works mainly 
focused on the Lyapunov exponents (LEs) and properties ob- 
tained from them (dimension, entropy, ...) |pl r7l— PlOn . Less is 
known about the associated tangent space directions, in par- 
ticular, there are no studies analyzing the structure of Lya- 
punov vectors (LVs) apart from the main one (pointing along 
the most unstable direction). Correlations of the main LV 
were only recently addressed iflll [l2ll taking advantage of the 
generic mapping of a DDS into an extended system, but with 
contradicting results. 

When considering LVs other than the main one, one must 
distinguish among different vector types. The so-called char- 
acteristic LVs 01311 constitute the only intrinsic set of vectors 
that is univocally defined and is covariant with the dynam- 
ics. However, their numerical computation is difficult in high- 
dimensional systems. We adapt to DDSs one of the methods 
to compute characteristic LVs [ 14]. 

In this paper we report on the generic properties that LVs 
exhibit in DDSs in the large delay limit. Our numerical and 
theoretical results indicate that the LVs (the main one and 
the others) behave in qualitative and quantitative terms like 
in one-dimensional dissipative systems with spatiotemporal 
chaos. 



II. TIME-DELAYED CHAOTIC SYSTEMS 



delay: 



dt 



(1) 



where y T = y(t—r). In a Lyapunov analysis we are interested 
in the evolution of infinitesimal perturbations Sy(t): 



dSy 
dt 



u5y + v5y T , 



(2) 



which is also a delayed equation with 5y T = 5y(t — r), and 
where u and v are functions: u(y,y T ) = d y J-, v(y,y T ) = 
dy T J-, If the DDS is chaotic an initial random perturbation 
becomes after some time aligned with the main LV (reaching 
a stationary state in a statistical sense), and the average expo- 
nential growth rate is the largest LE of the system Ai . Numer- 
ically, this is simply achieved by integrating Eqs. (Q]) and (0 
for long enough times. For non-leading LVs, more involved 
algorithms are needed (see below). 

We focus our study on systems in which delayed and non- 
delayed terms are separated: J-(y, y T ) = Q(y) + TI{Vt)- Our 
results are expected to be generic for this class of models. In 
our study, we have assumed the non-delayed part is linear 1 , 
Q(y) — —ay, and all nonlinearities appear in the retarded 
component 1Z. Many important time-delayed systems, in- 
cluding the Mackey-Glass (MG) JH] and Ikeda models ifTfjll . 
and optical delayed feedback systems 11711 . can be expressed 
in this mathematical form. In particular, we have studied 
numerically different nonlinear functions: the Mackey-Glass 
(MG), TZ(p) = bp/ '(1 + p 10 ), IS; and the nonlinear function 
TZ(p) = 6sin 2 (p — tp) that appears in some optical cryptosys- 
tems with delayed feedback lfl7l [l8ll . Almost identical results 
in qualitative and quantitative terms are obtained for these two 
systems. Therefore, for the sake of brevity, we choose to 
present the results only for the MG model, also used in ll ill . 

In our simulations we have integrated numerically Eq. (Q]) 
using a third-order Adams-Bashforth-Moulton predictor- 
corrector method 11911 . while linear equations, e.g. (0, have 
been integrated using an Euler method with the non-delayed 



DDSs may exhibit chaos even in the simplest situation in 
which the main variable is a scalar and its evolution is de- 1 We carried out some simu i ations with Q( y ) = - ay 3 / 2 that did not , 
termined by a delay-differential equation with one constant vealed relevant differences. 
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part integrated semi-implicitly. The parameters we used for 
the MG model are a = 0.1, b = 0.2 and the time step was 
dt = 0.2. In a numerical integration of a time-delayed sys- 
tem H, the temporal discretization makes the system finite- 
dimensional (with r/dt degrees of freedom) but this has no 
effect in the results whatsoever since an increase of time res- 
olution (i.e. decreasing dt) does not modify the largest LEs 
(and their associated vectors). 

In several previous works [l2ll it was found use- 

ful to map the DDS into an equivalent spatially extended dy- 
namical system of "size" t that evolves at discrete times 9 as 
follows. We express the continuous time as 

t = x + 9r (3) 

where x G [0, r) is the "spatial" position and 9 e Z is the dis- 
crete time. This spatial representation can be applied to both, 
the state of the system, y(t) — > y(x, 0), and the perturbations, 
5y{t) — > 5y(x, 9). The benefit of this spatial representation is 
that one can analyze the DDS with the tools available for spa- 
tially extended dynamical systems. We shall use the spatial 
representation of delayed systems here to analyze LVs and an- 
alyze the existence of long-range correlations as well as some 
other dynamical properties. 

III. CHARACTERISTIC VS. BACKWARD LYAPUNOV 
VECTORS 

There is some degree of confusion in the literature regard- 
ing the definition and computation of LVs. Many authors 
think of LVs as the orthonormal frame of vectors that results 
as a byproduct of computing the LEs with the algorithm of 
Benettin et al. H20fl . These vectors are usually called Gram- 
Schmidt or backward LVs |21] and they span the subspaces in 
tangent space that, at present time t, have grown at exponen- 
tial rates A n since the remote past. Backward LVs {b n (x, 9)} 
depend on the particular scalar product adopted for the Gram- 
Schmidt orthogonalization in Benettin's algorithm (though the 
subspaces they span are genuine). This and other undesired 
properties render the backward LVs unsuited to analyze prob- 
lems like extensivity or hyperbolicity questions. 

Ruelle and Eckmann 11131 12211 noticed time ago that one can 
define a different set of LVs, the so-called characteristic l2lll 
LVs {g n (x, 9)}, that are intrinsic to the dynamics and are uni- 
vocally defined (i.e. independent of how the scalar product 
is defined). Each of these vectors is covariant with the dy- 
namics, and the corresponding ?i-th LE is indeed recovered 
by either a forward or backward integration of an infinitesi- 
mal perturbation initially aligned with the n-th characteristic 
LV. Actually, the perturbation will be aligned with the LV at 
all times (hence the covariance). However for backward LVs 
(other than the main one, because b\ = gi) the correspond- 
ing LEs are only recovered integrating backward. Until very 
recently, no efficient numerical algorithms were available to 
compute characteristic LVs in large extended systems. 

In this paper we calculate the set of characteristic LVs by 
adapting to DDSs the method recently introduced by Wolfe 
and Samelson 11411 . In order to find the nth characteristic 



LV one has to compute the first n — 1 forward LVs in addi- 
tion to the first n backward LVs. Forward LVs {f n (x, 9)} are 
obtained like backward LVs but integrating the perturbations 
backward in time (from the remote future to the present time 
t). They have to be calculated using the transpose Jacobian 
matrix so that LEs are obtained with the usual ordering |[2lll : 
Ai > A2 > • ■ ■ . For our time-delayed system this procedure 
encompasses to integrate: 

_dSy^)_ =u§y (t) + v6y(t + T). (4) 
at 

The diagonal coefficient u = d w T{w, z)\ w=y{t) z=y{t _ T) is 
identical to that in Eq. (fJJ. The off-diagonal term has a curious 
structure stemming from transposing a Jacobian matrix with 
nonzero elements r temporal units below the main diagonal. 
Thus, in Eq. ©, we have i = d z T(w, z)\ w=y[tJ< _ T)z=y{t) . 
Notice that, as the integration runs backward in time (t — > 
—00), the minus sign in the left hand side is canceled. As 
occurs with backward LVs, periodic Gram-Schmidt orthonor- 
malizations are needed to avoid the collapse of all perturba- 
tions along the most unstable direction. Note also that, due to 
the delay, Eq. ([l} cannot be integrated backward, and hence it 
is necessary to store the trajectory y(t) at every time step in 
the computer to be used in the computation of Q. 

Once both backward LVs and forward LVs have been com- 
puted the characteristic LVs are easily calculated following 
the prescriptions in Ref. ltl4ll . One delicate point in delayed 
systems is that one must make sure that both forward and 
backward LV sets correspond exactly to the same time interval 
[t,t + r). 



IV. LONG-RANGE CORRELATIONS OF LYAPUNOV 
VECTORS 

In this section we study the form of both backward (Gram- 
Schmidt) and characteristic LVs corresponding to the largest 
LEs, their localization properties, and the existence of long- 
range correlations. 

LVs are strong localized objects, therefore, like in Refs. ITll 
[l2tl . it is very convenient to work with the associated fields 
obtained after a logarithmic transformation: 

h n (x, 9) = ln\g n (x,9)\ (5) 

(and likewise for backward LVs). We will refer to h n as the 
nth surface due to the similarities of its dynamics with that 
of the kinetic roughening of fractal surfaces in growth mod- 
els (see details below). Also recall that the norm is irrelevant 
because a LV only indicates a direction in tangent space, and 
this implies that the mean height of h n is irrelevant (only its 
profile fluctuations matter); this should be beard in mind for 
the theoretical analysis below. 

Figure Q2a) shows a snapshot of the characteristic LV 
surfaces hi, h%, and /13 corresponding to the three largest 
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FIG. 1: (Color online) (a) Snapshot of the surfaces corresponding 
to characteristic LVs n = 1,2,3 (r = 3277 t.u.). (b) Difference 
of 2nd and 3rd LV-surfaces with respect to the 1st LV-surface. The 
existence of plateaus evidences the piecewise copy structure of LV 
surfaces with respect to the 1st one. 

LEs 2 . We find that non-leading LV surfaces are approximately 
piecewise copies of the leading LV surface hi, see Fig. [32b). 
This 'replication property' is a highly nontrivial phenomenon 
that was originally discovered to occur in chaotic extended 
dissipative systems If23"l l24ll . We also emphasize that charac- 
teristic LVs exhibit a tendency to clusterize, contrary to back- 
ward LVs whose localization sites are scattered due to the im- 
posed orthogonality. Indeed, as exemplified by the snapshot 
shown in Fig. |TJa), the first three characteristic LVs localize 
(reach their largest magnitude) at the same point (completely 
overlapping on the leftmost region in this particular snapshot). 
This does not occur all the time but in an intermittent manner. 

Interestingly, these structural features- namely, replication 
and clustering- have recently been reported to occur generi- 
cally for LVs in chaotic spatially extended systems l23Tl25ll . 
This deepens in the analogy between DDSs and systems with 
extensive chaos in one dimension, which happens to hold even 
at the level of non-leading LVs. 

Regarding the existence of long-range correlations, we find 
that LV surfaces have a self-affine spatial structure at long 
scales, which translates into power-law correlations. A de- 
tailed analysis of the spatial structure can be best achieved by 
computing the Fourier transform of the surfaces at discrete 
times 9 given by h n (k,9) = -5= fj" exp(ikx) h n (x,9) dx, 

with wavenumbers k 6 [^r,J^], where r is the system 
size in this representation. For a given LV surface h n , 
the Fourier transform of the two-point correlator (h n (xo + 
x,6)h n (xo,9)) — (h n (xo,9)) 2 ,is the so-called structure fac- 
tor: 

S n (k) = (h n (k,0)h n (-k,9)) (6) 

where the brackets denote average over time 9 and realiza- 
tions. Therefore, the structure factor S n (k) directly informs 
about the nth LV surface correlations at scale 1/fc. 

In Fig. |2] it can be seen that the structure factor for the first 
LV-surface decays asymptotically as k~( 2a+1 \ with the so- 



2 For r = 3277 t.u. there are 168 positive LEs. The no-th LE vanishes, with 
n ~ 0.0514 r + 0.79. 



backward LVs 




characteristic LVs 




FIG. 2: Structure factors of LV-surfaces for the MG model with 
r = 3277 t.u. The curves from top to bottom correspond to 
n — 1,4,8,16,32. The insets show the dependence of the local 
slopes of different curves. For n = 1 the asymptotic slope is —2 
whereas for n > 1 the asymptotic slopes are ~ — 1 and ~ — 1.2 for 
backward and characteristic LVs, respectively. 



called roughness exponent a = 1/2 (in agreement with the 
result in 111 211 for the Ikeda model). Figure [2] also reveals the 
existence of a particular crossover wavenumber fc* « =^ for 
the nth LV. The exponents observed for k < fc* are — 1.2 and 
— 1 for characteristic and backward LVs, respectively. These 
numbers coincide with those previously reported in dissipative 
systems in one dimension ll23ll24ll . Our interpretation of these 
results is that for backward LVs the fc _1 dependence seems to 
be the consequence of residual correlations with a geometric 
origin in the orthogonality of the basis. The —1.2 exponent 
of characteristic LVs indicates they convey more information 
among distant parts of the system than backward LVs; which 
can be understood as a consequence of the covariance of char- 
acteristic LVs as they are consistent with both past and future 
evolutions [characteristic LVs with n > 1 are "saddle solu- 
tions" of (O forced by (Q]i, see [24]]. 
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V. MAIN LYAPUNOV VECTOR AND THE 
UNIVERSALITY CLASS QUESTION 

As we have seen in the previous section, the surface hi as- 
sociated with the main LV exhibits scale-invariant correlations 
in the k — > limit. This strongly suggests that this LV sur- 
face should belong to one of the universality classes of sur- 
face growth. The fact that the LV surfaces obey scaling laws 
and that systems with spatiotemporal chaos could be divided 
into a few universality classes, according to the scaling of the 
associated surfaces, has implications in our understanding of 
chaos in extended systems from a statistical physics point of 
view. Moreover, the generic replication and clustering of char- 
acteristic LVs along the main vector direction in extended sys- 
tems B23I - I25II (also observed here for time-delay systems, see 
previous section) makes it even more interesting to determine 
the universality class of the main LV, since this is expected 
to provide a great deal of information about the structure of 
space and time correlations of the LV corresponding to lead- 
ing as well as sub-leading unstable directions. 

For a wide class of one-dimensional spatio-temporal 
chaotic systems (see [12]) the statistical features of the main 
LV surface are well captured by the l<i-Kardar-Parisi-Zhang 
(KPZ) M27I1 stochastic surface growth equation: 

d t h(x, t) = X (x, t) + n(d x h) 2 + ud xx h, (7) 

where \ is a white noise. The KPZ equation defines itself 
an important universality class of surface growth. Regarding 
chaotic systems, in many of them the main LV surface belongs 
to the KPZ class. This includes coupled-map lattices, cou- 
pled symplectic maps, the complex Ginzburg-Landau model, 
Lorenz 96 model, among others 11121 l23l l24l 28]. Only LVs 
of anharmonic Hamiltonian lattices are known 12911 to exhibit 
correlations that are clearly inconsistent with KPZ exponents 
(a > a KPZ = 1/2). 

Statistical mechanics teaches us that, given the fact that 
KPZ represents a dynamical universality class, one would ex- 
pect that a large collection of different systems could belong to 
the KPZ universality despite their apparent differences. Only 
symmetries and conservation laws would determine the uni- 
versality class. That may explain why Hamiltonian (energy 
conserving) systems, in contrast to dissipative systems, do not 
generally belong to the KPZ universality. On this basis, DDSs 
were also proposed Hill to belong to a different universality 
class, namely the Zhang model, arguing that these systems 
break the x — > —x symmetry in the spatial representation. 
However, this conclusion was in contradiction with an earlier 
work 111211 where KPZ was postulated. Our aim in this section 
is to clarify this question by means of a theoretical analysis 
and extensive numerical simulations. 



A. Theoretical analysis: KPZ vs. Zhang equation 

Our starting point is the evolution equation in tangent space 
©, of which the main LV is the asymptotic solution. In the 
spatial representation (see Sec. [TI]i we perform the change of 



variables 9 = (t — x) /r, so that (O is now written as: 

d x Sy(x, 9)=u Sy{x, 9) + v 6y(x, 6-1), (8) 

where x G [0,r) is the 'spatial' position and 9 G Z is the 
discrete time. We have seen in the previous section that a de- 
scription in terms of surfaces instead of vectors themselves is 
more appropriate due to the strong localization of the latter. 
Moreover, to relate Eq. dS) with one of the stochastic partial 
differential equations modeling surface growth, we wish to 
approximate the discrete differences in 9 by a partial deriva- 
tive. In sum, we have to take two steps: 

(i) Transform to a surface: h(x, 9) = In \Sy(x, 9)\. 

(ii) Approximate 9 by a continuous variable: s(9) — s(9 — 
1) -)■ d e s. 

Note that these two steps do not commute. 



1. Linear Zhang model 

The treatment of Eq. (0 by Sanchez et. al. [11] proceeded 
with step (ii) before step (i). In more detail, after step (ii), 
Eq. <(8j becomes: 

d e 5y=-(l/v)d x 5y + (u/v + l)Sy, (9) 

and now transforming to the surface picture [step (i)] one has: 

d e h{x,6) = -(l/v)d x h + (u/v) + l, (10) 

so that we get an equation for the field h(x, 9). This equation 
was already analyzed in Ref. lITTll and we summarize some of 
its properties in the following. On the one hand, the random 
drift term (l/v)d x h gives rise to an effective diffusion term 
Dd xx h. On the other hand, the presence of v in the denomi- 
nator induces large fluctuations of those terms proportional to 
C = 1/v every time that v takes values close to zero. Usu- 
ally, the probability density function is algebraic at the lowest 
order: 

p(«)~K (H«i), (id 

with a > —1 to be normalizable (and a = in gen- 
eral). In turn, P(() is heavy tailed with large fluctuations 
P(C) ~ |e|~( 2+cr ) (this is indeed in agreement with data we 
have collected in our simulations of the MG model: a w 0). 
These reasoning led the authors of Ref. [11] to conclude that 
the main LV surface in DDSs behaves following ( TTOb and 
generically falls into the universality class of the linear Zhang 
model Q: 

d e h(x,e) = Dd xx h + c + e(M), (12) 

that describes surface growth driven by a heavy-tailed noise 
£(x,6) with a probability distribution P(£) ~ |£|-( x +m) for 
|£| > 1 and the index ^ = 2(1 + a) (> 2 if a > 0). 
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2. KPZ universality 

In this work we propose an alternative analysis of Eq. ([8]), 
where step (i) is carried out before step (ii). This would be 
more suited because one applies the nonlinear transformation 
(i) prior to the approximation (ii). Thus after defining the cor- 
responding LV surface, h(x,6) = \n\5y(x,9)\, Eq. (O be- 
comes 

d x h{x,9) = u±vexp[h(x,9-l)-h(x,9)], (13) 

where the choice of sign ± comes from the absolute value and 
is irrelevant for the arguments that follow. We approximate 
the difference h(x, 9 — 1) — h(x, 9) by the partial derivative 
—dgh. With no further approximations 3 and taking logarithms 
in Eq. ( fT3l we get: 

d e h(x,9) =hx\v\-ln\d x h-u\ (14) 

where the term r\ = \n\v\ is again a fluctuating noise-like 
source due to the chaotic character of the trajectories. In the 
neighborhood of v = the probability distribution of v is (fTTT >. 
and thus large values of \r]\ are exponentially rare: P(jj — > 
— oo) ~ e~^ <7+1 ''' n ', i.e. the noise is not heavy tailed. 

We can obtain a more intuitive equation by making use of 
the small gradient approximation, \d x h\ 1. Expanding (fT4l 
in the form: 

dgh(x, 9)= v - In \u\ + ^ + %^r- + 0[(d x h)% (15) 

where the prototypical quadratic term (d x h) 2 of KPZ appears. 
If u is fluctuating, the drift term (d x h)/u again yields (at large 
scales, after a spatial averaging) an effective diffusion term. 
In this case, Eq. ( [TBl l would lead to the KPZ Eq. (0. In con- 
trast, if u is a constant, like for instance in the case of the 
MG model, a diffusion term would also eventually appear far 
from the small gradient limit as a result of the absolute value 
in In \d x h — u\ in (fT4l . Finally, notice that the noise-like term 
ln(|u|/|u|), which is different from that in (I12t . will not lead 
to rare events and its distribution will be exponentially decay- 
ing in general. 

At variance with Sanchez et al. derivation ifTTll 
(cf. Sec. IV A p we have arrived here at the LV surface 
equation (fT5T l by expanding on the local slope d x h, which 
is generally expected to be a reasonable approach when the 
roughness exponent a < 1 since the spatial average over a 
window of extent £ should scale as (\d x h\)e ~ and goes 
to zero as we coarse-grain t 3> 1. 

B. Growth exponent 

In the field of growing interfaces it is customary to quantify 
the temporal features by looking at the growth of the surface 



3 Considering the time derivative dg h small enough so that the exponential 
in Eq. 1131 can be expanded and only the lowest order may be retained 
yields Eq. (ToJ. We propose here a different derivation that avoids this 
assumption and, remarkably, leads to a different result. 
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FIG. 3: (Color online) Evolution of the width from a random per- 
turbation for three different values of r. The lower panel shows the 
local derivative. For the largest value of the delay, r = 52429 t.u., 
the local derivative exhibits a plateau about 0.63, close to the theo- 
retical value for KPZ (2/3 = 2/3). The curves are averages over at 
least 1000 realizations. 



width at a given time when started from a flat initial perturba- 
tion [h(x, 0) ~ const.]. To do so we let random initial per- 
turbations 5y(x, 0) to evolve and measure the average surface 
width at time 9 as 



W 2 (e) = ^[h(x,6)-h(9)] 2 y (16) 

where h{9) = t _1 C h(x, 9)dx is the mean surface position 
and brackets denote averaging over different initial realiza- 
tions. The growth exponent j3 is defined as the exponent of 
the transient power-law growth: W 2 ~ 9 2 @, before the sur- 
face fluctuations saturate (for 9 <c 6* x ); while W 2 ~ t 2q , 
when saturation is reached (for 9 3> 6 X ) and the perturbation 
is virtually aligned with the LV: Sy oc gi =>• h = hi + const. 
The key point now is that f3 equals 1/3 for KPZ, while it is 1/4 
for the Zhang model with fi > 2 13111 . This makes the time ex- 
ponent j3 an excellent index to distinguish between KPZ and 
Zhang behavior (note that the spatial exponent a is the same 
in both models: a = 1/2). 
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FIG. 4: (Color online) Height-height correlation functions G q (l), 
q = 1, . . . , 5. Multiscaling is observed at small I. Above a certain 
characteristic size l c ~ 10 no multiscaling is observed. Note that l c 
is insensitive to t, and the curves overlap almost perfectly (except at 
large / due to finite-size effects). 



for all q. For truly extreme event dominated fluctuations, this 
length scale is expected to diverge (however slowly) with the 
system size, so that the fluctuation distribution is truly non- 
Gaussian in the thermodynamic limit. This slow divergence 
was indeed measured for models in the Zhang universality 
class I2r3l . In the case of time-delay systems, this would corre- 
spond to having I c (t) increasing with r. In contrast, we find 
that the characteristic length l c ~ 10 remains constant even 
after an increase of the delay of 16 times. Again, the sim- 
ulations by Sanchez et. al. ifllTl were carried out in systems 
with delays that were too short to obtain conclusive evidence 
on the dependence of l c with r. Certainly, we cannot rule out 
for sure an extremely weak dependence of l c on r, but if it 
exists it must be sub-logarithmic and well beyond the preci- 
sion that we can reach in our simulations. We conclude that 
multiscaling of the LV surface fluctuations in DDSs seems to 
be a short scale phenomenon that has no effect in the thermo- 
dynamic limit where the universality class is defined; in this 
case described by the KPZ equation. 



Figure [3] shows that W 2 grows with an exponent that, 
for large enough delays, progressively approaches the KPZ 
growth exponent /3k pz = 1/3, which strongly supports the 
KPZ asymptotic (r — > oo) scaling for chaotic DDSs. In 
the paper by Sanchez et al. Hill the largest time delay used 
t = 4000 t.u. resulted in an estimation j3 = 1/4. In the light 
of our simulations with much larger delays we conclude that 
the value of t used in Ref. Hill was insufficient to detect the 
true asymptotic scaling exponent. 

C. Multiscaling 

In addition to the growth exponent, Ref. Hill also invoked 
the existence of multiscaling of the main LV surface as an 
argument supporting its identification with the universality 
class of the linear Zhang model (fT2b . This model, contrary to 
the KPZ equation, leads to surfaces that exhibit multiscaling, 
i.e. strongly non-Gaussian tails, induced by extreme events. 
If we compute the gth-height-height correlation function for 
points separated a distance I, 

G q (l) = (\h 1 (x + l,9)^h 1 (x,e)\"y /q , (17) 

(where the overline denotes the average over x and the brack- 
ets denote a realizations average) one finds that G q (l) ~ l a ", 
and multiscaling exists if the roughness exponents a q depend 
on the index q. 

We have computed the qth order height-height correlation 
function G q (l), given by Eq. ( fTTl l, for several values of r 
(i.e. system sizes); see Fig. [4] The distribution of the sur- 
face fluctuations shows clear signs of multiscaling with local 
roughness exponents a q that depend on q for length scales be- 
low a typical scale of about l c ~ 10. These non-Gaussian fea- 
tures disappear at larger scales (I 3> l c ), where the a q « 1/2 



VI. CONCLUSIONS 

In this work we have implemented (to our knowledge for 
the first time) characteristic LVs in DDSs. Adaptation of the 
method proposed in [ 14] to this kind of systems — together 
with the computer capabilities nowadays available — allowed 
us to reach systems with fairly large delays, which serves to 
investigate the "thermodinamic limit" of these systems. Our 
results for the LVs coincide quantitatively with those obtained 
in extended dissipative systems [23, 24]. 

In addition we have revisited the question of which univer- 
sality class the main LV belongs to. After simulations with 
very large delays we may conclude that the main LV sur- 
face falls into the universality class of the Kardar-Parisi-Zhang 
equation. Our theoretical arguments support this conclusion 
as well. 

In sum, DDS are equivalent to extended dynamical systems 
in the sense that infinitesimal perturbations exhibit the same 
exponents characterizing spatiotemporal correlations. DDS 
have been traditionally considered to be different because of 
the lack of extensivity of the Lyapunov spectrum: the positive 
exponents approach zero as ~ 1/r [7] and the (Kolmogorov- 
Sinai) entropy saturates with r. However the identification of 
r with a size implies that comparisons should be done in tem- 
poral units of = t/r, and extensivity is then recovered. As 
X n t = X n T0 = A n 6, the redefined LEs A n = tA„ do not 
decay to zero with r, and the Lyapunov spectrum converges 
to a density in the thermodynamic limit. 
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